# Set of experimental functions

load.model = function(fname) {
    model <- read.table(fname)
    model <- colnames("id", "group", "age", "sex")
    model <- as.data.frame(model)
    model$sex <- factor(model$sex, levels=c("Female", "Male"))
    model <- model[order(model$id),]
    return(model)
}

test.batch.kendall = function(outputfn, pattern='', path='.') {
    library(bigmemory)
    library(ff)
    library(irr)
    library(Rniftilib)
    source("~/r-nitools/neuroimage.R")
    source("~/r-nitools/base.R")

    list.cors = test.read.rda(pattern, path)
    list.voxmaps = test.covert(list.cors)
    vec.ws = test.run.kendall(list.voxmaps)
    ni.write.vector(vec.ws, '../mask_8mm.nii.gz', outputfn)
}

ff.setdesc = function(ff.obj, desc.file) {
    desc = list(
        filename = filename(ff.obj),
        vmode = vmode(ff.obj),
        dim = dim(ff.obj),
        dimorder = dimorder(ff.obj),
        names = names(ff.obj),
        dimnames = dimnames(ff.obj)
    )
    
    dput(desc, desc.file)
    
    return(desc)
}

ff.getdesc = function(desc.file) {
    return(dget(desc.file))
}

ff.load = function(desc.file) {
    desc = ff.getdesc(desc.file)
    desc$update = F
    
    ff.obj = do.call("ff", desc)
    
    return(ff.obj)
}

r2ff.save = function(r.obj, file.prefix, ...) {
    ff.obj = as.ff(
        r.obj,
        filename=sprintf("%s.ff", file.prefix),
        ...
    )
    ff.setdesc(ff.obj, sprintf("%s.desc", file.prefix))
    return(ff.obj)
}

big.load = function(desc.file) {
    attach.big.matrix(dget(desc.file), dirname(desc.file))
}

big.rm = function(big.obj) {
    rm(big.obj)
}

test.loadffs = function(pattern='desc', path='.') {
    fns = list.files(path, pattern=pattern)
    num.fns = length(fns)
    
    list.ffs = vector("list", num.fns)
    
    for (i in 1:num.fns) {
        cat(sprintf("Loading element #%03i out of %i\n", i, num.fns))
        
        list.ffs[[i]] = ff.load(fns[i])
        close(list.ffs[[i]])
    }
    
    return(list.ffs)
}

# Read in a set of Rda files with correlation matrices holding the cormats as a list
test.rda2ff = function(pattern='Rda', path='.', ff.path='./ff') {
    # Get absolute path info
    cur.path = getwd()
    setwd(path)
    path = getwd()
    setwd(cur.path)
    setwd(ff.path)
    ff.path = getwd()
    setwd(cur.path)
    
    fns = list.files(path, pattern=pattern)
    split.fns = strsplit(fns, "[.]")
    num.fns = length(fns)
    
    list.cors = vector("list", num.fns)
    
    for (i in 1:num.fns) {
        cat(sprintf("Converting rda -> ff for subject #%03i out of %i\n", i, num.fns))
        
        tmp.path = file.path(ff.path, split.fns[[i]][1])
        
        load(fns[i])
        diag(cor.mat) = 1
        
        cor.ff = r2ff.save(cor.mat, tmp.path, overwrite=T)
        close(cor.ff)
        list.cors[[i]] = cor.ff
        
        rm(cor.mat)
    }
    
    return(list.cors)
}

# Read in a set of ff files with matrices and convert to a big matrix + save
test.ff2big = function(ff.path="ff", big.path="bigmemory") {
    ff.fns = list.files(ff.path, pattern="desc")
    base.fns = sub(".desc$", "", ff.fns)
    num.fns = length(ff.fns)
    
    list.cors = vector("list", num.fns)
    
    for (i in 1:num.fns) {
        cat(sprintf("Converting ff -> bigmemory for subject #%03i out of %i\n", i, num.fns))
        
        ff.obj = ff.load(file.path(ff.path, ff.fns[i]))
        open(ff.obj)
        
        list.cors[[i]] = list(
            desc=sprintf("%s/%s/%s.desc", big.path, base.fns[i], base.fns[i]),
            data=NA
        )
        
        dir.create(sprintf("%s/%s", big.path, base.fns[i]))
        tmp.mat = as.big.matrix(
            ff.obj[,],
            type = "double",
            separated = T,
            backingfile = sprintf("%s.bin", base.fns[i]),
            descriptorfile = sprintf("%s.desc", base.fns[i]),
            backingpath = file.path(big.path, base.fns[i])
        )
        rm(tmp.mat)
        
        close(ff.obj)
    }
    
    return(list.cors)
}

test.r2big = function(pattern='Rda', path='.', big.path='./bigmemory') {
    fns = list.files(path, pattern=pattern)
    num.fns = length(fns)
    
    prefix.fns = sub(".Rda$", "", fns)
    
    list.cors = vector("character", num.fns)
    
    for (i in 1:num.fns) {
        tmp.path = file.path(big.path, prefix.fns[i])
        tmp.desc = sprintf("%s/%s.desc", tmp.path, prefix.fns[i])
        
        if (file.exists(tmp.desc)) {
            cat(sprintf("Reading filename fir subject #%03i out of %i\n", i, num.fns))
            list.cors[i] = tmp.desc
        } else {
            cat(sprintf("Converting subject #%03i out of %i\n", i, num.fns))
        
            load(fns[i])
            diag(cor.mat) = 1

            dir.create(tmp.path, mode="770")
            tmp.big = as.big.matrix(
                cor.mat,
                type="double",
                separated=TRUE,
                backingfile=sprintf("%s.bin", prefix.fns[i]),
                descriptorfile=sprintf("%s.desc", prefix.fns),
                backingpath=tmp.path
            )
            rm(cor.mat)
            big.rm(tmp.big)
            
            list.cors[i] = tmp.desc
        }
    }
    
    return(list.cors)
}

ff.get.voxmap = function(v, list.cors) {
    num.subs = length(list.cors)
    num.voxs = nrow(list.cors[[1]])
    
    voxmap = matrix(NA, num.voxs, num.subs)
    for (s in 1:num.subs) {
        is.fopen = is.open(list.cors[[s]])
        if(!is.fopen) open(list.cors[[s]])
        
        voxmap[,s] = list.cors[[s]][,v]
        
        if(!is.fopen) close(list.cors[[s]])
    }
    
    return(voxmap)
}


# Takes list of subject correlation matrices and converts into a list of different voxels maps
test.submaps2voxmaps = function(list.cors, ff.path="./ff/voxmaps") {
    cat("-----\n")
    num.subs = length(list.cors)
    num.voxs = nrow(list.cors[[1]])
    
    list.voxmaps = vector("list", num.voxs)
    
    # Convert
    cat("Converting", num.voxs, "voxmaps\n")
    for (v in 1:num.voxs) {
        cat(sprintf("#%03i...", v))
        
        # Load
        if (file.exists(sprintf("%s/vox_%05i.desc", ff.path, v))) {
            ff.voxmap = ff.load(sprintf("%s/vox_%05i.desc", ff.path, v))
        }
        # Save
        else {
            # Get voxel maps across subjects
            r.voxmap = matrix(NA, num.voxs, num.subs)
            for (s in 1:num.subs) {
                is.fopen = is.open(list.cors[[s]])
                if(!is.fopen) open(list.cors[[s]])

                r.voxmap[,s] = list.cors[[s]][,v]

                if(!is.fopen) close(list.cors[[s]])
            }
            
            ff.voxmap = r2ff.save(r.voxmap, sprintf("%s/vox_%05i", ff.path, v), overwrite=F)
            close(ff.voxmap)
            rm(r.voxmap)
        }
        
        list.voxmaps[[v]] = ff.voxmap
    }
    
    cat("-----\n")
    return(list.voxmaps)
}

# Takes list of subject correlation matrices and converts into a list of different voxels maps
test.submaps2voxmaps = function(list.cors, big.path="./bigmemory/voxmaps") {
    cat("-----\n")
    num.subs = length(list.cors)
    tmp.mat = big.load(list.cors[1])
    num.voxs = nrow(tmp.mat)
    rm(tmp.mat)
    
    list.voxmaps = vector("character", num.voxs)
    
    # Convert
    cat("Converting", num.voxs, "voxmaps\n")
    for (v in 1:num.voxs) {
        cat(sprintf("#%03i...", v))
        
        base.fn = sprintf("vox_%05i", v)
        tmp.path = sprintf("%s/%s.desc", big.path, base.fn)
        
        # Load
        if (file.exists(tmp.path)) {
            ff.voxmap = tmp.path
        }
        # Save
        else {
            # Get voxel maps across subjects
            r.voxmap = matrix(NA, num.voxs, num.subs)
            for (s in 1:num.subs) {
                tmp.mat = big.load(list.cors[s])
                r.voxmap[,s] = tmp.mat[,v]
                rm(tmp.mat)
            }
            
            tmp.mat = as.big.matrix(
                r.voxmap,
                type = "double",
                backingfile = sprintf("%s.bin", base.fn),
                descriptorfile = sprintf("%s.desc", base.fn),
                backingpath = big.path
            )
            rm(tmp.mat)
            
            ff.voxmap = tmp.path
        }
        
        list.voxmaps[v] = ff.voxmap
    }
    
    cat("-----\n")
    return(list.voxmaps)
}

test.covert2 = function(list.cors) {
    cat("-----\n")
    num.subs = length(list.cors)
    num.voxs = nrow(list.cors[[1]])

    list.voxmaps = vector("list", num.voxs)
    
    # Create new cor mats
    cat("Creating new empty list\n")
    for (v in 1:num.voxs) {
        list.voxmaps[[v]] = big.matrix(num.voxs, num.subs, type="double")
    }
    
    # Convert
    cat("Converting", num.voxs, "voxmaps\n")
    for (v in 1:num.voxs) {
        cat(sprintf("#%03i...", v))
        for (s in 1:num.subs) {
            list.voxmaps[[v]][,s] = list.cors[[s]][,v]
        }
    }
    
    cat("-----\n")
    return(list.voxmaps)
}

# Runs kendall's W on a list of voxmaps
# Note:
# each element of the list is a matrix
# each matrix has columns with different voxmaps (e.g. from diff subject)
# each matrix has rows with different correlations for a given voxel
test.run.kendall = function(list.voxmaps) {
    library(irr)

    num.voxs = length(list.voxmaps)
    num.subs = ncol(list.voxmaps[[1]])
    
    vec.ws = vector("numeric", num.voxs)
    
    cat("Computing", num.voxs, "voxmaps\n")
    for (i in 1:num.voxs) {
        cat(sprintf("#%03i...", i))
        is.fopen = is.open(list.voxmaps[[i]])
        if(!is.fopen) open(list.voxmaps[[i]])
        
        vec.ws[i] = kendall(list.voxmaps[[i]][,])$value
        
        if(!is.fopen) close(list.voxmaps[[i]])
    }
    cat("\n")
    
    return(vec.ws)
}

# Degree
# Obtains the mean correlation between each voxel with the rest of the brain
test.degree = function(list.voxmaps) {
    num.voxs = length(list.voxmaps)
    num.subs = ncol(list.voxmaps[[1]])
    
    degree.mat = matrix(NA, num.subs, num.voxs)
    
    cat("Computing across", num.voxs, "voxmaps\n")
    for (i in 1:num.voxs) {
        cat(sprintf("#%03i..", i))
        is.fopen = is.open(list.voxmaps[[i]])
        if(!is.fopen) open(list.voxmaps[[i]])
        
        degree.mat[,i] = apply(list.voxmaps[[i]][,], 2, mean)
        
        if(!is.fopen) close(list.voxmaps[[i]])
    }
    cat("\n")
    
    return(degree.mat)
}

# Distance-based multiple regression
# Run adonis on the between-subject comparisons of each region/voxels connectivity maps
test.adonis = function(list.voxmaps, adonis.formula, adonis.model, ...) {
    library(vegan)
    
    adonis.formula = sprintf("maps ~ %s", adonis.formula)
    
    num.voxs = length(list.voxmaps)
    num.subs = ncol(list.voxmaps[[1]])
    
    list.adonis = vector("list", num.voxs)
    
    cat("Computing across", num.voxs, "voxmaps\n")
    for (i in 1:num.voxs) {
        cat(sprintf("#%03i..", i))
        is.fopen = is.open(list.voxmaps[[i]])
        if(!is.fopen) open(list.voxmaps[[i]])
        
        maps = as.dist(1 - cor(list.voxmaps[[i]][,]))
        list.adonis[[i]] = adonis(as.formula(adonis.formula), data=adonis.model, ...)
        
        if(!is.fopen) close(list.voxmaps[[i]])
    }
    
    cat("\n")
    
    return(list.adonis)
}

tmp = read.table("../test_model.txt")
test.model = data.frame(
    comt=as.factor(tmp$rs4680),
    drd4=factor(tmp$X7repeat, labels=c("no_7R", "has_7R")),
    ext=as.numeric(tmp$EXT),
    age=as.numeric(tmp$Scan.Age)
)

# Takes an adonis data structure
# and converts it into various vectors (actually matrices) of p-values
test.getadonis.pval = function(list.adonis) {
    num.elems = length(list.adonis)
    vec.effects = rownames(list.adonis[[1]]$aov.tab)
    num.effects = length(vec.effects) - 2
    vec.effects = vec.effects[1:num.effects]
    
    adonis.pvals = matrix(NA, num.elems, num.effects)
    names(dim(adonis.pvals)) = c("voxels", "effects")
    colnames(adonis.pvals) = vec.effects
    
    for (i in 1:num.elems) {
        for (j in 1:num.effects) {
            adonis.pvals[i,j] = list.adonis[[i]]$aov.tab$Pr[j]
        }
    }
    
    return(adonis.pvals)
}


ni.adonis.get.pval = function(vec.effects, maps.compare) {
    num.regions = length(maps.compare)
    num.effects = length(vec.effects)
    
    m.effects = matrix(NA, num.regions, num.effects)
    names(dim(m.effects)) = c("voxels", "effects")
    colnames(m.effects) = names(vec.effects)
    
    for (r in 1:num.regions) {
        for (e in 1:num.effects) {
            m.effects[r,e] = maps.compare[[r]]$aov.tab$Pr[vec.effects[e]]
        }
    }
    
    return(m.effects)
}
